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I. INTRODUCTION 

A domain is a region in a ferromagnetic material with the magnetization along a given direction. A magnetic 
material contains many domains with different magnetizations pointing in different directions in order to minimize 
the total magnetostatic energy. Regions with different orientations of their magnetization can be close to one another 
albeit with a boundary called a domain wall (containing typically about 10^ - 10^ atoms). 

Saturation occurs when all these regions align along some common direction imposed an external applied field, the 
total magnetization reaching its largest value Mg. 

The width of a domain wall is equal to tt^J A/K where A is the typical nearest neighbor Heisenberg exchange 
interaction and K the typical anisotropy constant (see Table QJ. Hence, a magnetic wall results from exchange and 
anisotropy, being thinner for higher anisotropy or smaller exchange (In Fe it is about 30 nanometers whereas in a 
hard material like Nd2Fei4B it is about 5 nanometers, only). Domain wall energy is given by A^/AK illustrating once 
again the competing role of exchange and anisotropy. 

For bulk materials, walls of the Bloch type occur whereas in thin films Neel type walls are encountered when the 
film thickness is close to the exchange length (defined by iex — \J AjK^ which is a few nanometers for ferromagnetic 
materials like Ni, Fe or Co, see Tabled). In the case of soft or amorphous materials characterised by a vanishing 
anisotropy constant one uses rather the magnetostatic exchange length defined by = \/~AJM^. In all cases, 
the wall width 5 is obtained from the exchange length via 5 = i:iex- 

A single parameter Q = iK/M"^ allows to discriminate between simple {Q < 1) and complex wall profiles {Q > 1) 
(see Malozemoff and Slonczewski 0). For example, in fig.QJa Bloch wall, belonging to the class {Q < 1) is depicted 
with the magnetization rotating in a vertical plane. 

Mathematically, a domain wall appears as a result of a non-linear two-point boundary value problem (TPBVP) 
since it separates two distinct regions with a well defined value of the magnetization. The TPBVP originates from a 
minimization of the total magnetic energy that contains in general a competition between the anisotropy and exchange 
energies. 

In this work, a general numerical approach based on the relaxation method is applied to the study of domain profiles 
in several geometries: bulk, wires and thin films. 

This report is organised as follows: In section 2, the numerical relaxation method is described; in section 3 we discuss 
Bloch walls, whereas radial walls in cylindrical wires are described in section 4. In section 5 Neel walls are described 
and finally section 6 contains a discussion and a conclusion. 

II. THE RELAXATION METHOD 

Traditionally, TPBVP are typically tackled with the shooting method. The shooting method typically progresses 
from one boundary point to another using, for instance, Runga-Kutta integration with a set of initial conditions 
attempting at reaching the end boundary. 

For regular Ordinary Differential Equations (ODE), simple shooting is enough to reach the solution. In more 
complicated ODE, one has to rely on double shooting also called shooting to a fitting point. The algorithm consists 
of shooting from both boundaries to a middle point (fitting point) where continuity of the solution and derivative are 
required. In certain cases, one even has to perform multiple shooting in order to converge toward the solution 
In the case of presence of singularities (within the domain or at the boundaries) the shooting method in all 
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its versions: simple, double or multiple does not usually converge. We find that it is the case also with do- 
main walls because of a rapid drop of the solution somewhere in the integration interval (due to the rapid 
change of the magnetization orientation in the wall). In this work, we develop, a new method to tackle the domain 
wall problem based on the relaxation method and find it quite suitable to handle relatively fast changes in the solution. 

The basic idea of the relaxation method is to convert the differential equation into a finite difference equation 
(FDE). When the problem involves a system of N coupled first-order ODE's represented by FDE's on a mesh of M 
points, a solution consists of values for N dependent functions given at each of the M mesh points, that is iV x M 
variables in all. The relaxation method determines the solution by starting with a guess and improving it, iteratively. 
The iteration scheme is very efficient since it is based on the multidimensional Newton's method (see Numerical 
recipes 0). The matrix equation that must be solved, takes a special, block diagonal form, that can be inverted far 
more economically both in time and storage than would be possible for a general matrix of size (MN) x (MN). The 
solution is based on error functions for the boundary conditions and the interior points. 
Given a set of N first-order ODE's depending on a single spatial variable x: 

^ = g^ix, yi, . . . , yN),j - 1, 2, ... iV (1) 
we approximate them by the algebraic set: 

= Ek = Vk - Vk-i - {xk - Xk-i)gk{xk,Xk-i,yk,Vk-i), fc = 2, 3, . . .M (2) 
over a set of M — 1 mesh points defining [xk-i, Xk] intervals with fc = 2, 3, . . . M . 

The FDE Ek provide N equations coupling 2N variables at the mesh points of indices fc — 1, fc. The FDE's provide 
a total of M{N — 1) equations for the MN unknowns. The remaining equations come from the boundary conditions 
i: 

At the first boundary xi we have: = Ei = B{xi,yi) 

At the second boundary X2, we have: = Em+i = C{xM,yM) 

The vectors Ei and B have ni non-zero components corresponding to the ni boundary conditions at xi. The vectors 
Em+i and C have n2 non-zero components corresponding to the n2 boundary conditions at X2, with ni -\- n2 = N 
the total number of ODE's. 

The main idea of the relaxation method is to begin with initial guesses of yj and relax them to the approximately 
true values by calculating the errors Ei to correct the value of yj iteratively. Relaxation might be viewed as a rotation 
of the initial vector (representing the solution) under the constraints defined by Ei. The evolution of the relaxation 
process, is obtained from solution-improving increments Ayt that can be evaluated from a first-order Taylor expansion 
of the error functions Ek- 

It is that expansion that results in the matrix equation possessing a special block diagonal form, allowing inversion 
economically in terms of time and storage (see ref. Q). 



III. BLOCH WALLS 



The energy of an uniaxial ferromagnetic material comprises anisotropy and exchange terms. An infinite volume is 
considered to exclude any shape related demagnetization energy. Exchange energy density is given by 0: 

2 dxi dxk 

where Einstein summation convention is used for repeated indices i^k,l — 1...3. The uniaxial anisotropy en- 
ergy is given by KijMiMj with i,j — 1...3. For simplicity, we assume a single uniform exchange constant A 
(see Table and a sole dependence on the x coordinate of all components of the magnetization M. We have 
M = {0, AIs smd{x), Ms cos 9{x)) (see fig.^. 0{x), the angle the magnetization makes with the z axis considered as 
the anisotropy axis. The sought profile is the function 9{x). Ms is the saturation magnetization when all individual 
magnetic moments in the material are aligned along the same direction. 

Integrating over all the volume, the total energy is given by: 

^ = J_j2^^) +yM,^}d. (4) 
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This can be rewritten as: 



[A{^) +Ksin^e]dx (5) 
2 ax 

The energy minimum is found by nulhng the variational derivative of E with respect to 6. We find: 



dx 



2 Csin6icos6i = (6) 



with C = . The Bloch wah profile is given by the solution to the above second-order ODE written as a system of 
two first-order equations: 



dyi 
dx 

dm 

dx 

where t/i = S{x) satisfies the boundary conditions: 



V2 (7) 
^ sin cos 2/1 (8) 



lim yi{x) — tt; lim yi{x) — (9) 

X — > — 00 X — >OG 

It is understood that the sharp transition between the 6 = {) phase and the 9 ^ tt \s behind the failure of all shooting 
methods. 

Using the relaxation method, we easily obtain the wall profiles for any value of ^ as displayed in fig. |5| 
Actually, the Bloch wall problem is single scale and we can without performing the calculation for every ^ value, do 
it for one value and then change the scale accordingly. This is done as follows: The width of the domain wall is given 
by 5 = I/VJ as explained previously. We perform a scaling transformation to the x coordinate as: x ^ x/5 turning 
the ODE into: 



J2q 

— ^ - sin6lcos6l = (10) 
dx'^ 

The exact analytical solution of the above equation given by: 6{x) = 2tan^^(e^^) is indistinguishable from the 
relaxation method results displayed in fig. |21 

Numerically, this means one can do the calculation for = 1 and later on rescale the x variable in order to get the 
solution for any value (arbitrarily large or small) of ^. Despite the power of the relaxation method, we noticed that 
when ^ < 10~^ or when ^ > 10^, convergence becomes difficult due to rounding and conditioning errors. 

That rescaling works for many types of walls except Neel wall where we have an additional scale controlling the 
profile (see section V). 



IV. RADIAL WALLS 



The energy density of an infinite cylindrical (see figEJ uniaxial ferromagnetic material comprising uniaxial anisotropy 
and exchange terms is given by: 



For simplicity, we assume a single uniform exchange constant A and a sole dependence on the radial coordinate r 
of all the magnetization components of M . Integrating over all a cylindrical volume of radius i?, the total energy is 
given by: 
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^-^iMTr^ +^]+i.cos^.}2..d. (12) 

where 9 is the angle the magnetization makes with the z axis (see fig.|2Il. a plays the role of a lattice parameter, the 
minimal core radius, regularising the integral (see for instance 0)- As in the Bloch wall case, we consider that 9 is 
a function of one spatial coordinate only (r in this case) . Since the anisotropy energy is given by: K cos^ 9 with K 
positive, the base plane (perpendicular to the z axis) is easy, meaning the minimum of anisotropy energy is obtained 
when 9 = n/2 (see figO)). 

The total energy minimum is found by nulling the variational derivative of E with respect to 9(r). We find: 



d^9 ^ld9 ^sm29^^ ^ ] - (13) 
dr^ r dr 2 

with ^ = ^ . The radial wall profile is given by the solution to the above second-order ODE (equivalent to system 
of two first-order ODE's like the Bloch case) with the boundary conditions: 



lim 9{r) = 0; lim 9(r) = 7r/2 (14) 

r — >a r — *i? 

The limits: a ^ 0; i? — > oo are taken afterwards. 

Using several values of f we obtain the radial wall profile in figEl Again, like in the Bloch case, there is a single 
length involved and it suffices in fact to solve the TPBVP for a single case ^ = 1 and rescale all variables accordingly. 
This is not the case of Neel walls as decribed in the next section. 



V. NEEL WALLS 



Neel realized that in a regime where the thickness of a ferromagnetic film becomes comparable to the Bloch wall 
width, a transition mode within the plane can lower the total energy decisively. Unlike the Bloch wall problem where 
only two energy components (exchange and anisotropy) exist balanced by a single length scale, the Neel wall problem 
incorporates two characteristic length scales. The new length arises from the competition with an additional energy 
component, the internal field energy. This has important physical, mathematical and numerical consequences. On 
the physical side, a very rich behaviour of Neel walls in thin films was shown recently in ref. @. 
Domain structures in thin inhomogeneous ferromagnetic films with smooth and small inhomogeneities in the exchange 
and anisotropy parameters yield very complex domain structures Domain walls are fixed near certain inhomo- 
geneities but do not repeat their spatial distribution. In addition there are metastable chaotic domain patterns in 
periodically inhomogeneous films. 

The mathematical description of Neel walls entails the introduction of an internal magnetic field H created by p the 
induced pole density induced by the rotation of M. Mathematically we have divH=-divM=p. The magnetization 
is expressed as M = [Ms sin6'(x), Mg cos 9{x),0) in the xyz coordinates defined in fig. |S1 Since the divergence of M 
is not zero, we have an induced pole density p. In contrast, p = in the Bloch wall case since we recall in this case 
(see section III), M = {0, MsSm9{x), Alg cos9{x)). Assuming as done previously, that the components of H depend 
solely on the spatial variable x (see fig. 0), we obtain: 

dM^ „ dism9) 

p = -divM = = -Ms^ — - (15) 

ox dx 

The ODE that controls the wall profile 9{x) is derived exactly as before (taking account of the exchange and 
anisotropy terms) with the addition of the Zeeman term accounting for the presence of the internal field H{x): 

2A—-Ksm29 + MsH(x)cos9 = (16) 
dx-' 

The uniaxial anisotropy term is K sin 29 with 9, the angle the magnetization makes with the y axis (the anisotropy 
axis) . 

Note that the demagnetization energy (due to the finite thickness of the film along the z direction) is zero, since it 
is given by 27riVijMjM, with iV^^, = Nyy = 0, N^^ = 1, and = 0. 
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The difference between this equation and the previous ones (Bloch and Radial cases) is that the internal field term 
H depends on the profile 6{x). Writing H{x) instead of H{9{x)) makes the wall-profile equations non- autonomous 
because of the explicit x dependence in H(x). Additionally these equations are integro-differential because of the 
dependence of H{x) on 9{x) (see for instance j^). 

In this work we consider the thin film approximation and retrieve a system of three ODE's by introducing a third 
fmiction = H{x)/Hk with Hk = 2K/Ms the anisotropy field. The ODE system to solve is written with respect 
to normalised coordinates x — x/5 where 5 is the wall thickness {5 ~ 1/%/^ where, as before, ^ = 



dyi 
dx 
dy2 
dx 
dm 
dx 



= y2 (17) 

= sin yi cos yi - ^3 cos j/i (18) 
= -7rC?/2 cos j/i (19) 



The magnitude of the coupling constant C — -^jf- has a strong effect on the solution of the system. In the limit 
C = we recover the simple case with no internal field H{x) = like the Bloch wall case. As the magnitude of C 
increases, we get a greater variation in the spatial dependence of ys^x) and the system might become unstable and 
display numerical oscillations in spite of a drastic reduction of the integration step. 

We convert the boundary conditions from the ] — oo, +00 [ interval to the [0, +oo[ interval: 

lim yi{x) 7r/2; lim yi{x) = lim ysix) = (20) 

X — '0 X — 'OO X — >oo 

We describe below a special algorithm, that we developed, based on the relaxation method coupled to an iterative 
procedure. The pseudo-code follows: 



1. Initially, we introduce a guess profile (say 0o{x)), extract from it the pole density poix) using ea.llSland determine 

dx 



from it the field derivative using the divergence equation: '^^J-^^ = p(x). 



2. The ODE system is solved and that allows us to extract a new profile (say 6*1 (x)) that yields a new pole density 
pi{x) (using ea. [TK|l . 

3. We repeat this procedure to the n-th step with a profile 9n{x) yielding a pole density pn+i{x) that provides 
a profile 6n-i-i{x). The procedure stops when the difference between the two profiles 9n{x) and 0„-|_i(a;) in the 
mean-square sense becomes smaller than an error criterion. 

The latter profile will have then relaxed self-consistently to the sought optimal profile that minimises the total 
energy (see ref. and references within). 

The results we obtain with various values of C for 9{x) and the internal field H{x)/Hk are displayed in fig. [T] and 
fig. El The analytical result obtained for C = (Bloch case), given by: 9{x) = 2tan~-'^(e~^) is displayed in fig. Inland 
is indistinguishable from the numerical results we obtain with the relaxation method on the svstem lT^ 
The results obtained for the internal field displayed in fig. [7| show, as expected (see for instance ref. Q), that when 
C increases, the field (absolute) amplitude becomes larger close to the origin. In addition, as C increases the field 
extends to larger distances farther from the origin. That, in fact, points to the origin of the integro-differential nature 
of the problem. Inspection of eq. ^| shows that in addition to the usual length scale (wall width) 6 = 1/ a/Cj we have 
another length given by: Sjsi = ^ A/ {KC) = 6/VC arising from the internal field whose strength is given by the 
coupling constant C. As C increases, non-local effects increase, the length ratio Sn/S = l/VC decreases (making the 
competition between the two lengths harder to deal with because of the disparity of the two lengths) and it becomes 
more and more difficult for the relaxation method to find an optimum result satisfying the TPBVP. 



VI. DISCUSSION AND CONCLUSION 



The magnetic domain profile is a challenging mathematical and numerical problem. In this work, we treated 
with the relaxation method, in the simple domain structure case (Q < I), wall configurations in several interesting 
physical cases: Bloch walls in ferromagnetic bulk systems, radial walls in cylindrical ferromagnetic wires and the Neel 
walls in thin ferromagnetic films assuming in all cases uniaxial anisotropy. 

In the Neel case, we showed than in the thin film approximation (in present technology, thin means ~ 10 — lOOA) one 
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is able to solve the wall problem with the relaxation method with a proper selection of the variables. Nevertheless, a 
major difficulty appears at higher value of the thickness t along the z direction (see fig. 13). 

When the thickness of the film increases the system becomes a full integro-differential system whereas in the thin 
film approximation, we get a set of coupled non-linear ODE's that we have to treat with a special self-consistent 
algorithm. The non-locality of the internal field is responsible for the appearance of logarithmic tails in the spatial 
variation of the magnetization angle. That means the TPBVP must be solved over an ever increasing interval size. 
The algorithm we have developed still applies but one has to use the finite thickness formulas for the field H{x) 

(eq. I24|l and its derivative ^^^j^ (eq. I25|l as shown in the Appendix. The extension of this work to other types of walls 
(originating from other types of anisotropy for instance, or the complex wall shape case Q > 1) or wall dynamics is 
challenging since the wall profile rapid change imposes a constraint on the time integration step. 
Previously, Smith treated domain wall dynamics in small patterned magnetic soft thin (~ lOO/im) films and turned 
the dynamic Landau-Lifshitz equations into a set of coupled non-linear ODE's. It turns out that the system of 
equations, he found is stiff (see for instance ref. imposing a very small integration timestep slowing down 

considerably the integration process on top of the difficulties of the TPBVP. 

The extension of this work to domain structures in inhomogeneous media (see ref. 0) is also quite interesting, 
particularly to the case of thin magnetic films that are of high technological interest such as recording, memories 
(Magnetic RAM's and Tunnel Junctions) and Quantum computing and communication. 
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APPENDIX 

We derive, in this Appendix, the formula for the internal field from the induced pole density. The magnetization 
is expressed as M = {AIsSm6{x),MsCOs9{x),0) in the xyz coordinates defined in fig. Since the divergence of M 
is not zero, we have an induced pole density p{x). The internal field is obtained from the pole density by integration 
accounting for the finite thickness of the film. 

Using divH — —p we infer from general theorems of elcctromagnctism that: 



H{r) = I dr'p{r')^ ^ (21) 

with r = {x, y, z),r' ^ (x' , y', z'). 
By symmetry we have Hy ~ Hz = and the x component H{x) in the plane z = is written as: 

Hix) = / dx' / dy' / dz'pix') — (22) 



47r 



*/2 [{x - xf + [y ~ yf + (zff 
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A first integration over y' gives: 



H{x) = -— dx' dz'pix') ^ , ' ^ (23) 

A second integration over z' yields the result: 

H{x) = - / p{x')t^n-\- -]dx' (24) 

ttJ^oo 2{x~x') 

The relation divJJ = p gives the integral expression of ^^^j^ needed in the integration of system of ODE's eq. ^1 

dH{x) 2t 



, p{x') 5 dx' (25) 

dx TT J-ac A{x~x'f+f^ 

In the finite thickness case, one needs the solve an integro-differential system of equations defined by the system 
of ODE's [T^ and the integral eq. In the case of thin films t ^ 0, we recover from eq. |2Slthe previous definition 
^^^Ip- = p{x) by using the 5 function definition: 

1 2t 

5(x -x') = - lim ^ (26) 

7rt^0 4(2;_2;')2+t2 

TABLES AND FIGURES 



Material 


T, 


fioMs 


A 


K 




Unit 


[K] 


[T] 


10-"[J/m] 


10^ [J/m^] 


[nm] 


Fe 


1044 


2.16 


1.5 


0.48 


2.8 


Co 


1398 


1.82 


1.5 


5 


3.4 


Ni 


627 


0.62 


1.5 


-0.057 


9.9 


Permalloy 


720 


1.0 


1.3 





5.7 


Cr02 


393 


0.5 


0.1 


0.22 


3.2 


SmCos 


993 


1.05 


2.4 


170 


7.4 



TABLE I: Properties of Ferromagnetic Materials: Tc is Curie temperature, /lo is vacuum permeability, Ms is saturation 
magnetization, A is exchange constant, K is magneto-crystalline anisotropy constant and iex is exchange length. Note that in 
the case of Permalloy (Ni^iFeioo-i: alloys with x ~ 80), one uses the magnetostatic exchange length defined as lex = \/ A/M} 
since A" ~ 0. 



FIG. 1: Behaviour of the magnetization direction for a Bloch wall. For an arbitrary point along the x axis, the magnetization 
M whose rotation is entirely confined within the vertical zOy plane makes the angle 6 with the vertical z axis, the anisotropy 
a:xis. 




FIG. 2: Variation of the magnetization angle with distance for a Bloch wall for various values of the exchange anisotropy ratio 
^. The analytical result 6{x) = 2tan~^(e~'^), for ^ = 1, is indistinguishable from the relaxation method result. 
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FIG. 3: Cylindrical coordinates displaying the spatial variation of the magnetization angle with radial distance from the wire 
axis. The base plane perpendicular to the wire axis z is an easy plane. 
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FIG. 4: Variation of the magnetization angle with radial distance r from the wire axis for various values of the exchange 
anisotropy ratio ^. As we increase f the angle increases faster from (Easy axis along z) to 7r/2 (Easy plane A-z) 
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FIG. 6: Variation of the magnetization angle with distance for a Neel wall for various values of the coupling constant C. 
Uppermost curve is for C = 0, whereas lower curves correspond respectively to C=l, 10, 50 and finally 100. The exact result 
corresponding to zero thickness along the z direction is indistinguishable from the C = curve. 
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FIG. 7; Variation of the normalised internal field H{x)/Hk with normalised distance for a Neel wall for various values of the 
coupling constant C. Uppermost curve is for C = 0, whereas lower curves correspond respectively to C=l, 10, 50 and finally 
100. 



